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Abstract: 

In this paper, we consider performance analysis of the decentralized power method for the eigendecomposition of the sam¬ 
ple covariance matrix based on the averaging consensus protocol. An analytical expression of the second order statistics of the 
eigenvectors obtained from the decentralized power method which is required for computing the mean square error (MSE) of 
subspace-based estimators is presented. We show that the decentralized power method is not an asymptotically consistent estima¬ 
tor of the eigenvectors of the true measurement covariance matrix unless the averaging consensus protocol is carried out over an 
infinitely large number of iterations. Moreover, we introduce the decentralized ESPRIT algorithm which yields fully decentralized 
direction-of-arrival (DOA) estimates. Based on the performance analysis of the decentralized power method, we derive an analyti¬ 
cal expression of the MSE of DOA estimators using the decentralized ESPRIT algorithm. The validity of our asymptotic results is 
demonstrated by simulations. 

keywords: decentralized eigendecomposition, power method, decentralized DOA estimation, ESPRIT, averaging consensus. 


1. Introduction 


Centralized processing in sensor networks requires the collection of measurements or sufficient statistics from all 
sensor nodes at a fusion center (FC) before processing to obtain meaningful estimates. A major drawback of a such 
centralized processing scheme with a single FC, is the existence of communication bottlenecks in large sensor networks 
with multi-hop communications fi~4}[20) . Averaging consensus (AC) protocols |3f7j|8j j24[|25| achieve an iterative fully 
decentralized calculation of the average of scalars that are distributed over a network of nodes. AC protocols use only 
local communications between neighboring nodes, thus, avoiding multi-hop communication. Moreover, AC protocols 
perform computations at the nodes and require no FC. Thus, AC protocols eliminate communication bottlenecks. 
These attributes of AC protocols make them attractive and a fully scalable alternative to centralized processing schemes 
in large sensor networks G3- 

The eigendecomposition of the sample covariance matrix is required in many applications, such as signal detection 
121 -23), machine learning |TJ, and DOA estimation [|9 TT] 13] 13} TS). Conventionally, the eigendecomposition is 
carried out in a centralized fashion, which hinders its application in large sensor networks. In fl4| , an algorithm which 
achieves a fully decentralized eigendecomposition of the sample covariance matrix is introduced. This algorithm 
combines the conventional power method (PM) |5} p. 450], which represents a centralized iterative eigendecomposition 
algorithm, with the AC protocol to achieve a fully decentralized eigendecomposition of the sample covariance matrix. 
We refer to this algorithm as the decentralized power method (d-PM). Analytical expressions of the second order 
statistics of the eigenvectors and eigenvalues of the conventional (centralized) sample covariance matrix are presented 
in [2} Theorem 9.2.2], The expressions from |2j| are asymptotic in the number of samples, i.e., they become accurate 
as the number of samples increases. In [28] a different approach, which is asymptotic in the effective signal-to-noise 
ratio (SNR), is proposed. This approach holds even for the case of one sample if the noise variance is sufficiently 
small and can be used to derive a general performance bound for DOA estimation 129 30j|. However, the accuracy 
of the eigendecomposition obtained from the sample covariance matrix using the d-PM not only suffers from finite 
sample effects and finite PM iterations, but also depends on the convergence speed of the AC protocol. This additional 
mismatch is introduced by the decentralized implementation and can be mitigated if the AC protocol is carried out 
over a large number of iterations. However, a large number of AC iterations is not always possible since it is associated 
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with a large communication overhead and latency. A performance analysis of the decentralized eigendecomposition 
which considers estimation errors introduced by the AC protocol is of wide interest for a large variety of applications. 

In we presented a fully decentralized DOA estimation algorithm using partly calibrated arrays. Our DOA 
estimation algorithm combines the d-PM with the conventional ESPRIT algorithm [ [13| and is thus referred to as the 
decentralized ESPRIT (d-ESPRIT) algorithm. The numerical simulations carried out in [ |20| and {19) show that the 
d-ESPRIT algorithm achieves similar performance as the conventional ESPRIT algorithm when a large number of AC 
iterations is used. However, an analytical study of the performance of the d-ESPRIT algorithm which supports these 
simulations has not been considered so far. Moreover, the behavior of the d-ESPRIT algorithm when only a small 
number of AC iterations is carried out has not been studied before. 

The first and main contribution of this paper consists in the derivation of an analytical expression of the second 
order statistics of the eigenvectors for the sample covariance matrix computed using the d-PM. Based on this expres¬ 
sion, we show that the d-PM is not a consistent estimator of the eigenvectors of the true measurement covariance 
matrix, unless the AC protocol is carried out over an infinitely large number of iterations. Moreover, we show that 
when the number of AC iterations used in the d-PM converges to infinity, our expression and the conventional expres¬ 
sion in j2] Theorem 9.2.2] become equivalent. The second contribution of this paper consists in the derivation of an 
analytical expression of the MSE in DOA estimation using the d-ESPRIT algorithm. 

The remainder of this paper is organized as follows. The measurement model is introduced in Section [2] In 
SectionRl we briefly revise the AC protocol f24) and the d-PM | fl4] and their main properties, used later in the analysis. 
Section B] considers the performance analysis of the d-PM, namely the computation of the second order statistics of 
the eigenvectors resulting from the d-PM. The d-ESPRIT algorithm and its performance analysis are presented in 
Section[5] Simulation results in Section[6]compare the performance of the d-PM and the d-ESPRIT algorithm with our 
analytical expressions of Sec.[4]and Sec7|5]and with the Cramer-Rao Bound (CRB) for partly calibrated arrays ©■ 

In this paper, (-) T , (•)* and (■) n denote the transpose, complex conjugate and the Hermitian transpose operations. 
The element-wise (Hadamard) product, diagonal or block diagonal matrices and the trace of a matrix are denoted by 
0 , diag[-] and Tr[-], respectively. The real part and the argument of complex numbers are denoted by Re[-] and arg[-], 
and j is the imaginary unit. The expectation operator and the Kronecker delta are expressed as E[ ] and v while 
[-]i, 1 1, 0, and lj stands for the (*,j)th entry of a matrix, the ith entry of a vector, the ixi identity matrix, the all 
zeros vector of size i and the all ones vector of size i, respectively. The sets of real and complex numbers are denoted 
by R and C. Variable x at the r/th PM iteration is expressed as X'‘ l> . The conventional (centralized) estimate of x is 
denoted by x whereas its decentralized estimate (computed using the AC protocol) is denoted by x. In decentralized 
estimation, we distinguish between estimates which are available at all nodes, such an estimate at the kth node is 
denoted by xua, and estimates which are distributed among all the nodes, where the kth node maintains only a part of 
the corresponding estimate denoted by x k . 

2. Measurement Model 


We consider a network of M = sensors clustered in K nodes, where the kth node comprises of M k sensors. 

The assignment of the sensors to the individual nodes is characterized by the sensor selection matrix T, whose entries 
are defined as 


J 1, if the ith sensor belongs to the y th node 
0, otherwise, 


( 1 ) 


where i = 1 ,,M and j = 1,..., K. 

The measurement vector at the kth node at time t is denoted as z k (t) £ C Mkx 1 and the overall measurement 
vector is denoted as 


z(t) = [z 1 {t),...,Z K (t)\ £ 


-<Mx 1 


( 2 ) 


The random measurement vector z(t) is assumed to be zero-mean with covariance matrix R = E [z(t)z H (t)]. The 
eigendecomposition of the true measurement covariance matrix R is defined as 


R^UAU h , 


(3) 


where U = [tii,... ,Um] and A = diag[Ai,..., Am] and U \.... ,Um are the eigenvectors of R, corresponding to the 
eigenvalues Aj > ■ • ■ > Am- For the later use, we define the matrix 

Bi±U-iTi 1 U h _ u (4) 

where[70; = ... ,u M ] and Ti = diag[Ai - A.., Xi-i - A;, A z+ i - A;,..., A M - A;]. 
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In practice, the covariance matrix R is not available, and can only be estimated from N observations of z(t),t = 

1.. .. ,N as 

1 N 

R=—^ Z {t) zH (t). (5) 

V t= 1 

We refer to the conventional estimator of the sample covariance matrix in Eq. M as the centralized estimator, since it 
requires that all measurements from every node are available at a FC. Let A, U, A, and u, be the estimates of A, U. A, 
and u, for i = 1 ,..., M, respectively, obtained from the centralized eigendecomposition of the sample covariance 
matrix R. 

In the following section, the decentralized estimation of the eigenvalues and eigenvectors of the true covariance 
matrix R using AC protocol, i.e., without a FC, is introduced and analyzed. We denote as A, E7, A, and u, the decen¬ 
tralized estimates of A, U, A* and Ui for i = 1..... Af respectively. 

3. Averaging Consensus and the Decentralized Power Method 

In this section, the AC protocol and its convergence properties are reviewed. Moreover, the decentralized eigendecom¬ 
position using the d-PM fj~4) is revised. 

3.1. Averaging Consensus 

Let X \,..., Xk denote K scalars which are available at I\ distinct nodes in the network, where the fcth 

the fcth scalar. Denote the conventional average of these scalars as x = A x k- I' 1 AC protocols 

is computed iteratively, where at the /Ah AC iteration, the fcth node sends its current local estimate of the average x'j!' L> 

to its neighboring nodes, denoted as the set A4> And receives the corresponding average estimates of the respective 

neighboring nodes. Then, the fcth node updates its local estimate of the average as follows 

(p) {p- 1) , (p- 1) /cx 

X fc = Wk,kX)e + 2^ W i,k x i (6) 

i&Mk 

where Wi^ is the weighting factor associated with the communication link between node i and node fc, which satisfies 

Wi t k = 0 when i £ A4 ^24], The AC iteration in Eq. ([(ijl is initialized with = Xk for fc = 1,..., K. For more 
details, see |24|. 

Denote with W the matrix whose entries are [W],; j = Wij for i, j = 1,..., K and let x^ = [x^f\ ..., x^] T , 
then, the update iteration in Eq. ([6]) can be expressed as 

x (p) = Wx &- x) = W 2 x (p ~ 2) = ■■■ = W p x w . (7) 

Iteration (j7]i converges asymptotically (for p —> oo ) to the vector of averages xIk if and only if 

W p (8) 

Let the eigendecomposition of the matrix W be 

W\p 1 ,...,P K ] = \fii, ■ ■ ■ ,Pk\ diag[ai,..., a K \, (9) 

where ot\ > • • • > o k. According to [24|, the matrix W which satisfies the asymptotic convergence condition <[8]) 
possesses the following properties: 

PI: The principle eigenvalue of the matrix W is unique (single multiplicity) and equals to one, i.e., o ( = 1. The 
corresponding normalized principal eigenvector of the matrix W is given by /31 Jk 1k - 

P2: The remaining eigenvalues of W are strictly less than o-\ in magnitude. 

In the following, we assume that the weighting matrix W satisfies the convergence condition ([8ji, which permits the 
use of properties PI and P2 in our analysis in Sec. [4] We express the decentralized estimate of the average x at the fcth 
node using p AC iterations as 

i [k ] = [W p x(%, (10) 

where W p x !{r> ]k denotes the fcth entry of the vector W p x^. The notation xm is used since the corresponding average 
is computed using the AC protocol and every node stores locally the computed average. 


node stores only 
7.8, 24 II 25 I, x 
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3.2. The Decentralized Power Method 


In this section, we first review the conventional (centralized) PM [5J p. 450], then we review the d-PM fl4) . 

The conventional PM is an iterative algorithm which can be used to compute the eigendecomposition of the sample 
covariance matrix R. Let us assume that l — 1 eigenvectors of R have been computed using the PM, then the Zth 
eigenvector is computed using the iteration 


_ 


H 


= ( Im - Ui-iUi- ± )Ru 


;(?-!) 


(ID 


where u q) denotes the Zth eigenvector of R at the r/th PM iteration, Im is the M x M identity matrix and Ui-\ = 
[ill,..., is the concatenation of the Z — 1 previously computed eigenvectors of R. The vector iq°' is a random 
initial value. If the PM is carried out for a sufficiently large number of PM iterations Q, then, the normalized vector 


Ui=U 


[ Q) /II«[ Q) | 


( 12 ) 


is the Zth eigenvector of the matrix R. 

The d-PM |jl4|] performs the computations in Eq. ( fTT] > and Eq. ( p~2| ) in a fully decentralized fashion based on the 

, where the kth node 


AC protocol. The key idea of the d-PM is to partition the Zth vector as u)' 11 = [u^ T , 

■ (4) ^ X 1 , 7 .-,(9) 


l l,K 


stores and updates only the fcth part, tij £ C" tX1 , of the vector uj HJ . Note that the notation Ui± is used, since the 
vector Ui t k is computed using the AC protocol and stored only at the fcth node. In the d-PM, iteration (111 is split into 
two steps. In the first step, the intermediate vector 


,; {a) ± ru\ q - 1] 


(13) 


is calculated. In the second step, the vector u\ q ^ is updated as 


■M) 


".'(9) fr tt H ,-/(9) 
u l— lU l— 1^1 5 


(14) 


where [7)_i = Ui- 1 ], In the following, we review how both Steps (13 i and (14 1 and the normalization Step 

( [T2l > can be carried out in a fully decentralized fashion 
Substituting Eq. Q into Eq. ( [13] ), yields 


u 


'(q) 


N 




(9) 
l ’ 


(15) 


where z t q l > =z H (t)u l <] l> . Note that =J2k-i z k(^)^fk w h ere z k is computed and stored locally at the 

kth node. Thus, in analogy to (JToj), the estimate of zfi at the fcth node computed using the AC protocol is 


~(q) _ is 

z t,i,[k] ~ 


W p [z?(t)u\ q - 1 \...,z»(t)u < i q - 1) ] 


- k 


(16) 


where P is the number of AC iterations used in this protocol |:14|. Using N parallel instances of the AC protocol, the 
kth node will locally maintain the scalars Thus, each node k can locally compute one part of the vector 


~f(q) ~/(g) 

ur ) as 


= jt z k(t) z ti rfci’ ^at in turn perform the first step of the d-PM iteration described in Eq. (| 13). 


Note that in the second step of' the d-PM iteration only the second term of Eq. (fl4|) has to be computea in 


a 


decentralized fashion 14 


H 


. This term can be written as where u- q> = In 

analogy to (|15|l, each node can locally compute its corresponding part of f//_ i (/(_ , U/ 1 ' 1 ’ once the scalars {u- 'j* }-Z-] 


are available at every node. This can be achieved using Z — 1 parallel instances of the AC protocol as 


<?[*]= * 


W Pl [uf, u {q) 


f.H f/(9)l T 
i, 1 u l, 1 > • ■ ■ i u i,K u l,K 


J k 


(17) 


where j is the zth scalar computed at the fcth node and Pi is the number of AC iterations used in these l — l AC 
protocol instances. Thus achieving the second step of the d-PM iteration. 
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After a sufficiently large number of PM iterations Q, the vector u) 
tion can be carried out locally once the norm 
as 


(Q) 


is normalized as in Eq. (12 1 . This normaliza- 


IK (Q) H[ fc ]=^ 


is available at each node which is achieved using the AC protocol 

(18) 


w p 2 [««>*««) 


*u 


-(Q)H 
H,K 1 


,(Q)it 

1,K 


J k 


where P 2 is the number of iterations used in the AC protocol instance. Thus, using the d-PM, the eigendecomposition 
of the sample covariance matrix can be calculated without FC. 

In the d-PM, communication between the nodes is required to compute the scalars in ( fi~6| ), ( | 1 7[ ) and (18 1 . From 
a signaling perspective, the first and most expensive computation is that of the N scalars in Fq. (] 16[>, where N AC 
protocol instances, i.e., PN AC iterations, are carried out to compute these scalars. The second most expensive 
computation lies in Eq. (17 1 which requires l — 1 AC protocol instances. The third and least expensive computation is 
the normalization of the eigenvectors which requires only one AC protocol instance. 


4. Performance Analysis of The Decentralized Power Method 

In this section, we first reformulate the d-PM as an equivalent centralized PM. Based on the centralized formulation, we 
derive an asymptotic analytical expression of the second order statistics of the eigenvectors for the sample covariance 
matrix obtained from the d-PM. Moreover, we show that the d-PM is not a consistent estimator of the eigenvectors of 
the true measurement covariance matrix R. 


4.1. Assumptions 


Our performance analysis focuses on the errors resulting from using a finite number of AC iterations P < oo to 
compute the scalars in Eq. (16 1 , because, from a signaling perspective, this step represents the most 

expensive calculation in the d-PM. Thus, tne following assumptions are made: 


A1: The number of AC iterations Pi and P 2 used to compute the scalars in (17 i and the normalization factors in (18 1 , 


respectively, are large compared to the number of AC iterations used to compute the scalars iff/ 1 K=i ’ '- e ^ 
Pi ©> P and P 2 P. Thus, errors resulting from the finite number of AC iterations in Eq. ( | 18 [ > and Eq. 0 
are negligible compared to those in Eq. ( fj~6| ). 

A2: The number of PM iterations Q is sufficiently large such that the errors resulting from the finite number of PM 
iterations are negligible. 


4.2. Error Expressions for the Decentralized Power Method 


The decentralized eigendecomposition of the sample covariance matrix using the d-PM yields the vectors 
Since under Assumptions A1 and A2 these vectors depend on P and not on I\ , P 2 and Q, we denote them as 
{ui{P)}i£ 1 - Due to finite AC iteration effects (P < oo), these vectors do not exactly correspond to the eigenvectors 
of the matrix R. The following theorem provides further insights into the properties of the vectors {'«/ (P) } . 

Theorem 1. Under Assumption Al, the vectors {ui(P)}ff. 1 are the eigenvectors of the matrix 

R(P) = K ( TW p T t ) © R , (19) 


where T is the sensor selection matrix defined in Eq. (|7| and R is the centralized sample covariance matrix defined in 
Eq. (|5f. 

Proof. We prove Theorem[l]by induction. Thus, first we prove that the vector U \ (P), which is computed using the 
d-PM, is the principle eigenvector of the matrix R{P). Then, assuming that the vectors {u fiP)} 1 ^ are the principle 
l — l eigenvectors of the matrix R(P), we prove that iii (P) is the Zth eigenvector of the matrix R{P). For convenience, 
we drop the dependency on P from R(P) and ufiP), throughout the derivations. 

Note that when the d-PM is used to compute the vector u\, then Eq. ([T4| reduces to 


iZ (9) - - 
Ul ~N 


N 


;(?) 

ATM 


zT(t), 




i T 


’ Z ip,[K] Z K (t) 


( 20 ) 
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since the matrix Uq = 0. Let 


Z(t) = 


(zi(t) 0 Ml 

0 M 2 Z 2 {t) 


\0 m k 0 m k 

where Zk(t) is defined in Eq. (J2ji. Then, Eq. (20) is written as 


UMi 

°Mk_ 1 

ZK{t)) 


-.(g) 





u 


(9-1) 


( 21 ) 




(9-1) 


t=l fc=l 


where the eigendecomposition of the matrix VP in Eq. (|9ji is used, T is the sensor selection matrix defined in Eq. ([T]) 
and P is the number of AC iterations used to compute the scalars {zfl r fe i}£Li- The product Z(t)Pk can be rewritten 


as 


Z{t)P k = [Pk,iZi ( t ),..., Pk,KZK{t )\ J 
= (Tp k ) © z{t ) 


( 22 ) 


where 0 k = [/3fc,i, • • •, 0k, k] ■ Substituting Eq. (22 1 into Eq. (21 1 , yields 

N K 


-M 


-( ((m)o*(t)) (( Tp k )®z(t)) H 


,N 


u 


(9-1) 


= ( 4EE a ^ ( r PkPkT T )e,(z(t)z H (t)) ) w 


£=1 fc=l 


AT 


( 9 - 1 ) 

i 


£=1 fc=l 
K 


= (kJ 2^(T0kPZT T )Q±- J2 (z(t)z H (t)) 


-,(9-1) 


fc=1 


=(k(tw p t t )qr 


N ; 

-,(9-1) 


Thus, the decentralized computation of it! using the d-PM can be written as the following iteration 


u 


(9) _ 




(9-1) 


(23) 


where i? is defined in Eq. (19 1 . Note that Eq. (23 1 corresponds to the update procedure of the conventional PM applied 
to the matrix R. Thus, after a sufficiently large number of PM iterations Q, the resulting vector u\ (p converges 
(if normalized) to the principle eigenvector of the matrix R. It follows from Assumption A1 that the decentralized 
normalization of ii ^ is accurate. Thus, under Assumption Al, the vector resulting from applying the d-PM to the 
sample covariance matrix R is the principle eigenvector of the matrix R computed using the conventional PM. This 
concludes the first part of the induction. 

For the second part of the induction, we assume that the vectors {u;}*Z! are computed using the d-PM and they 
are the first l — 1 eigenvectors of the matrix R. Then, we prove the induction for the vector U/. 

The computation of the vector ill using the d-PM is achieved as follows. First, the vector which is defined in 
Eq. (13 1 is computed in a decentralized fashion. In analogy to Eq. (20 1 , can be rewritten as 

-,(9-1) 


-'(g) 


Ru, 


(24) 


Second the scalars {'©'E, = { are computed in a decentralized fashion. Since under Assumption Al the AC errors 
resulting from this computation are negligible, the decentralized iteration used to compute the vector it; is reduced to 


u 


(9) 


H 


= ( I M - Ui-JI^Ru) 


-,(9-1) 


(25) 
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Note that Eq. ( |25j ) is equivalent to the conventional PM iteration ( fTT) applied to compute the /th eigenvector of the 
matrix R. After Q iterations of the d-PM, the resulting vector uf' !> is normalized. Again under Assumption Al, the 
normalization is accurate, thus, the resulting normalized vector iii is the 1th eigenvector of the matrix R computed 
using the conventional PM. □ 

Theorem [T] shows that, when the d-PM is used with a finite number of samples N and a finite number of AC 
iterations P to estimate the eigenvectors of the true covariance matrix R, then three different types of errors 

occur: 

El: Errors resulting from the finite number of AC iterations P. These errors are expressed in the matrix (TW p T t ) . 

E2: Errors resulting from the finite number of samples N. These errors are expressed in R. 

E3: Errors resulting from the finite number of PM iterations Q , which we neglect as stated in Assumption A2. 

Since the averaging matrix W is assumed to satisfy the convergence condition ([8b, we conclude that KTW p T t —> 
ImIm as P —» oo, consequently, R(P) —► R. We remark that R(P) —► K (TW^T t ) © R when N —» oo, i.e., for 
a finite number of AC iterations P, the eigendecomposition of the sample covariance matrix using the d-PM is not an 
asymptotically consistent estimator of the eigenvectors of the true measurement covariance matrix R. 

Theorem |T] simplifies the performance analysis of the d-PM, since it provides a link to an equivalent centralized 
algorithm formulation which can be analyzed using the conventional statistical analysis techniques and results |2j|. In 
the sequel, we start our performance analysis by introducing the error vectors which represent El and E2 types of 
errors. Then, we compute analytical expressions for these errors and finally we derive the second order statistics of 
the eigenvectors obtained from the d-PM. 

For the centralized eigendecomposition, the sample estimate of the /th eigenvector iii of the true covariance matrix 
R is expressed as 

ui = ui + dui, (26) 

where the error vector Sui accounts only for the finite sample effects, i.e., E2 type of errors, used in the computation 
of the sample covariance matrix. The decentralized estimate of the /th eigenvector is expressed as 

u^^m + SiniP), (27) 

where the error vector SufiP) accounts for errors resulting from the finite number of samples and the finite number of 
AC iterations, i.e., El and E2 type of errors. Similarly, we define 

R = R + 6R (28) 

and 

R(P) = R + 6R(P), (29) 

where SR accounts only for E2 type of errors and SR(P) accounts for both El and E2 types of errors. Using the 
aforementioned notation, the second order statistics of the eigenvectors computed using the d-PM are expressed as 

E [Sui(P) 6u*(P)] and E[<5fi,(P) Su T m (P)}. 

The following theorem simplifies the computation of E[5iq(P) 5u^ l (P)] and E[<5£q(P) Su^ n (P)] by introducing a 
simple expression for Sui(P). 

Theorem 2. Under Assumptions Al and A2, the error vector Sui (P) is given by 

Sui (P) = Bi (SRm + hi (P)) 


where Bi is defined in Eq. 0. 

K 

hi (P) =KJ2^ T k Rt"ui , (30) 

k =2 

T k = diag[T7?fc], T is defined in Eq. (| 7 |) and f3 k and a k are defined in Eq. ([p]). 

Proof. See Appendix |A| O 
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Note that in Theorem [2] the El type of errors are expressed in terms of the vector hfiP) which depends on the 
number of AC iterations P and on the eigenvalues and eigenvectors of the weighting matrix W, except for the principle 


eigenvalue and eigenvector. Since the magnitude of ap- is strictly less than one for k = 2,..., K (see Sec. 3.1 1 , it 


follows from Eq. <30 1 that hfiP) —> 0 as P —> oo, i.e., the AC protocol is carried out for an infinitely large number 
of iterations. Consequently, SufiP) contains no El type of errors when P —t oo. In Theorem [2] the E2 errors are 
expressed in terms of the matrix SR. If an infinite number of sample is available, i.e., N —> oo then SR —» 0. 
Consequently, Sui(P) —> 0 as both P and N tend to infinity. 

Based on Theorem^ analytical expressions for E[<$ii;(P) Su^(P)] and E [SufiP) Su^ n (P)} are introduced in the 
following theorem. These expressions are useful for computing the MSE of estimators which are based on the d-PM 
as we show later for the d-ESPRIT algorithm. 


Theorem 3. Under Assumptions A1 and A2 


M 


E [Sii(P)S<(P)]=^J2 K 


UiUp Si 


^ ^ (At - A*) 

i^l 

BMP)h”{P)B* 


and 


where Si 


E [Su l (P)Sul(P)\ = 


A/Ay, 


uiui. 


- 1) 


N (A*-A m ) 2 

+ BMP)hl{P)Bl 

is the Kronecker delta, N is the number of samples, Bi is defined in Eq. i 


• and hi ( P) is defined in Eq. ( 30 1. 


Proof. See Appendix [B] □ 

Note that only the second terms in E[6«;(P) 5u^ l (P)] and E[5«;(P) Su^fiP)] depend on the number of AC it¬ 
erations P (through the vectors hfiP) and h m (P)) and as P —> oo, these terms converge to zero. Consequently, 
E[(Ju;(P) SUm(P)} and E[<5iq(P) SuJ n (P)] tend to the centralized case found in j2j, when P — > oo. Moreover, as 
N —► oo for P < oo, E[<5iq(P) Su^fiP)] and E[f>{q(P) Su^ n {P)] do not converge to zero, i.e. the d-PM is not a 
consistent estimator for unless P is infinitely large. Theorem [5] shows that, in the second order statistics 

of the eigenvector estimates, the AC errors appear as an additive error term which is remarkable since in Theorem[l] 
the corresponding errors for the sample covariance are expressed as an element-wise multiplication with the matrix 
( TW p T t ). 

Note that in practice P can not be chosen to be arbitrarily large, thus, the second terms in E[£«;(P) Su^ l (P)} 
and E[dTi/(P) SuJ n (P)] will always be non-zero. However, P can usually be chosen such that the second terms in 
E[5ti;(P) Sullf P)] and E[<5iq(P) 5u^ n (P)\ are of the same order as the first terms. The proper choice of P will be 
further addressed in the simulations in Sec. [6] 

5. Performance Analysis of the d-ESPRIT Algorithm 

In this section, we briefly review the decentralized ESPRIT (d-ESPRIT) algorithm presented in |20) . Then, results 
from Theorem [3] are applied to derive an analytical expression for the MSE of DOA estimation using the d-ESPRIT 
algorithm. 

5.1. Signal Model and the ESPRIT Algorithm 

Consider a planar sensor array composed of I\ identically oriented uniform linear subarrays, where the Pth subarray 
is composed of Mp sensors. The distance d between two successive sensors measured in half wavelength is identical 
for all subarrays, see Fig. |T| The displacements between the subarrays are considered to be unknown, thus, the array 
is partly calibrated. 

Signals of L far-field and narrow-band sources impinge on the system of the subarrays from directions 6 = 
[0i,..., 6lY ■ The output of the fcth subarray at time t is given by 

Zk(t) = A k (0)s{t) +«fe(f), (31) 
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where Zk(t) = [z k \ (t)...., z k ^M k {t)} T is the baseband output of the M k sensors, s(t) = [si(i),..., sl(T)] t is the 
signal vector of the L Gaussian sources and n k (t) = [n* i(t),..., n kt M k ( t)] T is the vector of temporally and spatially 
complex circular white Gaussian sensor noise. The steering matrix of the fcth subarray is A k {6) = [a k (0\),... ,a k {6i)\, 
where a k (()i) is the response of the fcth subarray corresponding to a source from direction 9i relative to the array broad¬ 
side 

a fc (0,) = e^ K ‘ 

where m = [sin 0;, cos Oi} 1 and £ k is the position of the first sensor in the fcth subarray relative to the reference sensor 
in the first subarray, which is considered to be unknown. The measurement vector of the array is 

z(t) = A(8)s(t) + n(f), (33) 

where z(t) = [zf (t),... ,z^(t)] T , as in Eq. <[ 2 }, A{0) = [A\ (0),... ,Aj(0)] r and n(f) = [nf (t),... ,nj < (t)] T . 
The measurement covariance matrix is 


1 jird sin 0; 

X , c . 


„j(M k -l)Trdsin6i 


(32) 


R = E[z(t)z H (f)] = A(6)PA H (0) + a 2 I M (34) 

where P = E[s(f)s H (i)] is the source covariance matrix and a 2 is the noise variance and E[n(t)n H (t)] = <j 2 Im■ The 
eigendecomposition of the measurement covariance matrix can be partitioned as 


R = U S A S U? +U n A n U” 


(35) 


where A s £ R LxL and A n £ l )*( m l ) are diagonal matrices containing the signal and noise eigenvalues, 

respectively, U s = [«i,... ,Ul] £C MxL and U n = [ul+ i, ■ • • ,«m] £C Mx ( M_i ) are the signal and noise eigenvector 
matrices, respectively, and U\,... . u m are the eigenvectors of the matrix R corresponding to the eigenvalues Ai > 
• ■ ->Al> Al+i = . • • = Am=ct 2 , see Eq. 


The ESPRIT algorithm exploits the translational invariance structure of the measurement setup. This invariance 
structure is expressed in Fig. |T] where the sensors of the array are partitioned into upper and lower groups. The upper 
group contains the first M k — 1 sensors of the fcth subarray and the lower group contains the last M k — 1 sensors of 
the fcth subarray. The two groups are related by a translation with lag d 1131. Let J = diag[Ji,..., Jk] and J = 
diag[J-!,... ,J K \, where J k — [IM k ■ OmJ and J = [0M k -i,lM k -i], denote the selection matrices corresponding 
to the upper and lower groups, respectively. Based on the selection matrices, we define two matrices U s = JU S and 
t/ s = JU S . In the conventional Least Squares ESPRIT 
matrix 

W = 


13), the DOAs are computed from the eigenvalues of the 


(ufu.) 


-1 _ 


ufu . 


(36) 


as follows 

Oi = sin -1 (arg ) / (nd )) (37) 


where ^ for l = 1,..., L are the eigenvalues of the matrix \t, see [131 for details. 

In practice, the true covariance matrix |34| > is not available and its finite sample estimate R is calculated from N 
snapshot of the array output as in Eq. (Jsjl. Let U s , f/_,. C/ s , u,. Tf and ft be the estimates of (7 S , U s . U„-u,. ’P, and ■)/);, 
respectively, obtained from the eigendecomposition of the sample covariance matrix R. 


5.2. The d-ESPRIT Algorithm 


The decentralized ESPRIT (d-ESPRIT) algorithm which is presented in 1201 comprises two decentralized steps, first, 
the decentralized signal subspace estimation using the d-PM, second, the decentralized estimation of the matrix . 

The decentralized signal subspace estimation is carried out as explained in Sec. 3.2 The resulting decentralized 
estimate of U s , denoted as U s , is distributed among the subarrays, where each subarray stores only a part from Z7 S . 
In [201, based on the AC protocol, a decentralized algorithm for estimating the matrix W is introduced. Denote the 
corresponding estimate at the fcth subarray as ^[k]- In [20 , the computation of is achieved by rewriting Eq. (36i 
as 

= F[k]> (38) 
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'Y y'Y 
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Figure 1: Array topology and sensor grouping for 6 subarrays. The read lines indicate communication links between 
neighboring subarrays. 


where C\ fc i and Fw] are respectively the decentralized estimate of ufu s and uflJ 3 , at the fcth subarray. The AC 
protocol is used to compute each entry of the matrices C\j~] and F^-:\ as follows: 


,{k] 


=K 


W P3 [u'\ J | J, u' r K J K J ku :] ,k ]T 


- k 


(39) 


and 


fi.,j, [k]—K W'^U^ J 1 , . . . ,U^ K J K J K Uj^K} J 


(40) 


where and denote the (i,j)th entries of the matrices C^} and f^], respectively, and P 3 denotes the 

number of AC iterations used to compute Yj.[fc] and Thus, the fcth subarray can estimate the matrix as 

presented in Eq. (38 1 . The DOA estimation is carried out locally in the fcth subarray using the eigenvalues of the 
matrix ®[k] in Eq.lm 

To simplify the performance analysis of the d-ESPRIT algorithm, the following assumption is made 

A3: The number of AC iterations 1\ which is used to compute 'P / . from U s is large compared to the number of AC 
iterations used to compute U s , i.e., P 3 P. 


Under Assumption A3, the AC errors in the decentralized estimate of 'F are negligible compared to those in U s . Thus, 
the decentralized estimate of d/ is a good approximation of the centralized one using U s , i.e., 


® [fc] = ® = I u a u 


H 


■ H . 


u s u s , 


(41) 


where U s = JU S , U s = JU S and fc = I..... A. Let 'tl’i tor l = 1,..., L be the eigenvalues of (E'. In the d-ESPRIT 
algorithm, i/ji is used as an estimate of wi. Thus, we define the estimation error S'tPi as 

ipl-i’l + Hu (42) 

for l = 1..... fc, where S'tPi accounts for both El and E2 types of error. 


5.3. The MSE for DOA Estimation based on d-ESPRIT 

In p2| , the MSE of DOA estimation using the conventional Least Squares ESPRIT is presented. Assumption 
A3 allows us to use the analysis from (12 by replacing and Efdu; <Ju^] with K[Siii(P) Su^P)] and 

E[<5£q(P) Su^ n (P)], respectively. Thus, the following result from [ 121 holds. 


E[(<50 ; ) 2 ] = 


E[|^| 2 ]^Re[(^) 2 E[(^) 2 ]] 
2 (irdcos 9i) 2 


(43) 


where 


h , 


M6iln\ 2 ]='Yf l E[SU a rrf t SU a fr z 


(44) 
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( 45 ) 


E[(S'ip i y}=^E[6U s nr{6U s ]#*?, 


U S =U S + SU s , 7 P 4 q, ( Ufu s ) t/f (J - yH) and /xf 4 (t/f Ef s ) t/f (J - ^J). The It h left and 
righ t eig envectors which correspond to the /th eigenvalue of the matrix \D are denoted as qi and ff, respectively. In 
Eq. (44 1 and Eq. (45 1 , we replaced the conventional error SU s by SU s in the corresponding expressions of 121. 

Using the expression from Theorem[3] the expectation of the right hand side of Eq. ([44]) is rewritten as 


¥.[6U S rirf Suf] =EEK] i /[^( p ) a f( p )] 

XiXj 

(A, - A,-) 2 


L L 


2=1 j = 1 
L M 


= V EE 


H 


i=l i=1 

*=1 j=i 

Similarly, the expectation of the right hand side of Eq. (j45j is written as 

L L 

E[^s nrf 5U T s ) =Y.Y. [ r ' r H ntulP) 6uj (P)} 

2=1 j = l 
L L 


^ [r*/f j XiXjUiUj 


3 

3 4 * 


NiXi-XjY 


i —1 3=1 


(46) 


(47) 


Equations (|43]>—(|47|) provide the analytical expressions of the MSE for the DOA estimator using the d-ESPRIT 
algorithm. The second terms of Equations ( |46| > and ( [47] ) differ from the expressions in (I2| . Note that because of 
these terms, the MSE does not approach zero even if an infinitely large number of samples is available, i.e. the d- 
ESPRIT algorithm is not a consistent estimator of the DOAs. However, the simulations in Sec. [6] demonstrate that for 
a finite number of samples and a moderate SNR, a finite number of AC iterations is sufficient to achieve a performance 
comparable to that of the conventional ESPRIT algorithm 03- and to achieve the CRB [ |T6| . 


6. Simulation Results 

An array composed of AT = 6 subarrays each containing 2 sensors, i.e., M = 12, separated by half-wavelength is used 
in the simulations. The locations of the first sensors at the subarrays are (0,0), (0.45,0.99), (3.02, 0.45), (5.61, 0.90), 
(8.03,1.46) and (8.70,0.50) measured in half-wavelength. The upper and lower selection matrices of the 7 th subarray 
are J/, = [1,0] and J k = [0, 1 ]. The array topology depicted in Fig. [T] is assumed. Thus, the neighboring sets are 
N\ = { 2 ,3}, A /2 = {1,3}, A /3 = { 1 , 2 ,4}, A /4 = {3, 5, 6 }, A /5 = {4, 6 } and A/e = {4, 5}, where the fcth subarray 
communicates only with its neighbors A4- The entries of the weighting matrix W are selected as follows 

f max{card [A7 ] ,card \Mj ]} ’ ^ * Y 3 and J ^ M 

W\i-J = Ui, if * = J ( 48 ) 

{(). otherwise, 

where card[A/i] is the number of elements in the set A/}. The weighting factors {wi\f =1 are chosen as Wi = 1 — 
J2f= 1 w i,j’ re f er to 1241 for further details. This choice of the weighting factors only requires that each node knows 
the degree of its neighbors, thus, local but not global knowledge about the network topology is required at the node 
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level. The weighting matrix W resulting from the weighting scheme in Eq. (48 1 guarantees asymptotic convergence 


of the AC protocol, provided that the graph associated with the network is not bipartite 

Signals from L = 3 equal-powered Gaussian sources impinge onto the array from directions —14°, —10° and 5°. 
In the sequel, we evaluate our analytical expressions for the performance of the d-PM and the d-ESPRIT algorithm. 


6.1. Performance of the d-PM 


In the first set of simulations, shown in Fig. [2] and Fig. [3] the performance of the d-PM is evaluated as follows. 
We estimate L = 3 eigenvectors of the sample covariance matrix at the ith realization using the d-PM, i.e., we 
compute U s (i) = [«i(z),U 2 (*), 113 (f)], for i = 1 ,..., 200 realizations. Since the eigendecomposition is unique up 
to a multiplication with a unity-magnitude complex scalar, we use the method introduced in |4j Eq. (54)] to compute 
this scalar and correct the estimated eigenvectors. Then, the error matrix 6U s (i) = U s (i) — U s is computed, where 
U s = \u\ is the true signal subspace. At the vth realization, we define the normalized square error (SE) of the 

d-PM as 

SE d .pM(i) = Tr[6U s (i)6Uf (i)]/Tv[U s U”]. (49) 

Finally, we compute the root mean square error (RMSE) 


200 


RMSE d .p M 4 


200 


^ SE d _ PM (i) 


(50) 


The RMSE which is obtained from our analytical expression for the d-PM algorithm is denoted as ARMSE d . PM and 
it is computed as 

ARMSE d . PM 4 (^E[<5ti,(P)<5^(P)]/TV[17 s 17f] ] , (51) 


1 


where E[<5iq(P) <5u,^(P)] is given in Theorem^ 

In Fig. [2] we compare RMSE d .pM from 200 realizations with the ARMSE d .pM at different SNRs where the number 
of samples is fixed to N = 100 and the number of AC iterations P is taken to be 10, 20 and 30. The number of PM 
iterations is fixed to Q = 10 for all simulations. It can be seen from Fig.[2]that the error in the estimated eigenvectors 
decreases with increasing SNR until it reaches a certain value, which depends on P, then it is saturated. Note that 
the error computed using the analytical expressions ARMSE d _pM corresponds well to the one computed over 200 
realizations RMSE d .pM- 

In Fig. [3] the SNR is set to 10 dB and RMSE d .pM and ARMSE d _pM are computed for different numbers of samples 
N for three different numbers of AC iterations 10,20 and 30. The number of PM iterations is fixed to Q = 10. From 
Fig. [3] it can be observed that the error in the estimated eigenvectors decreases with N for small values of N. However 
when N is large, RMSE d _pM and ARMSE d _pM do not change with N as it can be seen in Fig. [3] for P = 10 and 
P = 20. For P = 30, RMSE d _pM and ARMSE d _pM show a good correspondence at very large values of N (which is 
not displayed in the figure). Moreover, a larger number of AC iterations results in a smaller error. This behaviour of the 
RMSE d _pM is in accordance with our conclusion that the d-PM is a consistent estimator of the eigenvectors of the true 
measurement covariance matrix only when P is infinitely large, see Sec. [4] It can also be observed in Fig. [3] that the 
error computed using the analytical expressions ARMSE d _pM corresponds to the one computed from 200 realizations 


of RMSE d .pM- 


6.2. Performance of the d-ESPRIT Algorithm 

In the second set of simulations, whose results are shown in Fig. [4] and Fig. [5] the performance of the d-ESPRIT 
algorithm is evaluated and compared to the analytical expressions of Sec. [5] In these simulations, we assume that the 
number of sources L = 3 is known to all subarrays, which is the case in many applications, e.g. communications, 
condition monitoring and acoustics. In applications where L is not known eigenvalue-based detection criteria available 
in the literature (such as the MDF |22j and the approaches in j27) and |26]) can be adapted in the decentralized scenario 
to detect the number of sources. The RSME of the d-ESPRIT algorithm is computed over 200 realizations as 

RMSE d . ESPRIT 4 ( 200 E 3 £$« ^ ’ (52) 
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Figure 2: The performance of the d-PM as a function of SNR for a fixed number of samples N = 100. 


where 9i(i) is the estimate of 9i computed at the /th realization using the d-ESPRIT algorithm. The analytical expres¬ 
sion of the RMSE of the d-ESPRIT algorithm is 


ARMSEd.ESPRiT — 


, i=i 


E[(<%) 5 


(53) 


where E[(56*;) 2 ] is computed using Equations (43 i—(47 i. In this set of simulations, the number of PM iterations is set 
to Q = 2. Moreover, in this simulation we plot the RMSE of the conventional ESPRIT algorithm computed as in 
Eq. (|52|, along with the corresponding performance analysis of the conventional ESPRIT algorithm from fl2) , which 
we denote as ARMSEesprits and the CRB for the conventional partly calibrated arrays 113- 

Fig. [4] demonstrates ARMSE^-esprit and RMSE^esprit for different SNRs where a fixed number of samples 
N = 100 is assumed. Note that at low SNRs the performance of the d-ESPRIT algorithm is similar to that of the 
conventional ESPRIT algorithm and it improves with increasing SNR. However, at high SNRs, it can be observed 
that the performance of the d-ESPRIT algorithm deviates from that of the conventional ESPRIT algorithm. It is clear 
from Fig. |4] that this deviation depends on the number of AC iterations P. Thus, for P = 30 and SNR values up to 
SNR = 15 dB the performance of the d-ESPRIT algorithm is similar to that of the conventional ESPRIT algorithm 
and both achieve the conventional CRB, whereas for P = 10 this deviation starts from SNR = 0 dB. Moreover, it 
can be seen from Fig. [4]that the RMSE of the d-ESPRIT algorithm at high SNRs is saturated and cannot be decreased 
unless the number of AC iterations is increased. 

In Fig. [5] the SNR is fixed to 10 dB and ARMSEj-esprit and RMSEj-esprit are computed for different number 
of samples N. It is obvious that the error in the d-ESPRIT algorithm does not approach zero when N —> oo, which is 
in accordance with our conclusion in Sec. [5] that the d-ESPRIT algorithm is not a consistent estimator of the DOAS, 
unless the number of AC iterations P is infinitely large. 

In Fig. [4] and Fig. [5] it can be observed that the values obtained for the averaged RMSE of d-ESPRIT algorithm 
HMSHd-KSPRiT are similar to the results of the analytical expression ARMSE^esprit- 


7. Conclusions 

In this paper, we derived an analytical expression for the second order statistics of the eigenvectors for the sample co- 
variance matrix computed using the d-PM. This analytical expression is used to derive the MSE of the DOA estimates 
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Figure 3: The performance of the d-PM as a function of N for a fixed SNR = 10 dB. 


obtained from the d-ESPRIT algorithm. It has been shown that the AC errors in d-PM and d-ESPRIT algorithm are 
dominant when N is very large or the SNR is very high. In our analysis, errors resulting from a small number of PM 
iterations is not considered. Nevertheless, in the simulations, it has been shown that 10 PM iterations are sufficient to 
make these errors negligibly small. 


A. Proof of Theorem [ 2 ] 


3.1 


Proof. In order to prove Theorem [ 2 ] the matrix R(P) is written in terms of R and W. Then, a first order analysis is 
carried out. For convenience, we drop the dependency on P from R{P), ui{P) and hfP) throughout the proof. 

The largest eigenvalue of the matrix W is a\ = 1 and its corresponding eigenvector is /3 1 = -^=1 k, see Sec. 
thus 

R = K ( TW p T t ) © R 

= K (rJ2a p k p k ^T T ^ QR 

= R + I< (E akTPkPZT^ ©R (54) 

K 

= R + Kj2^ k f k RT 

k—2 




where T k = dia.g[T/3/ ;; ] and, for the last equality, the rank one Hadamard product property j6| p. 104] is used. Note 
that the second term in Eq. ([54]) accounts for the errors resulting from the finite number of AC iterations P < 00 , and 
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Figure 4: The performance of the d-ESPRIT algorithm as a function of SNR for a fixed number of samples N = 100. 


that R —»• R when P —> oo. Substituting Eq. (28 1 in Eq. (54 1 yields 

K 

R = R + SR + K^2 a kT k (R + 5R)t 


k—2 

K 


(55) 


R + SR + K^ o£t k RT\ 


k—2 


where the term Yhk =2 °k Tk SRT k is neglected in the approximation since the magnitudes of 012 , ■ ■ ■, otK are all 
strictly smaller than one (see Sec. 3.11 and they are multiplied with the small variation SR. 

Multiplying Eq. ( |55j ) from the right with Ui and keeping the first order terms, we find 

Riii « |R + SR + K ^ a k T k RT k ^ (m + Su) 

« Rui + SRui + hi + RSii , 


■H 


(56) 


where hi is defined in Eq. (30 1 . The left hand side of Eq. (56 1 can be written as Rui = XiUi, where A/ is the hh 
eigenvalue of R. Expressing A/ as a perturbation in A/, i.e.. A; = A; + SXi, the left hand side of Eq. (56 1 becomes 

Riii = (a 1 + SXi^j {ui + Siii) 

~ A m + SX1U1 + A iSiii, 

where only first order terms are kept. Substituting Eq. ( |57j ) in Eq. ( |56| ) yields 

(R - XJ M ) Siii ~ SXiUi - SRui - h t 
The matrix R — XiIm can be written as 


(57) 


(58) 


M 


R — XiI m = ^ (Afc — A i)uiu k 

k=l,k^l 

= U-iT-iU H l 


H 


(59) 
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Figure 5: The performance of the d-ESPRIT algorithm as a function of N for a fixed SNR = 10 dB. 


where Z7_/ and r_/ are defined in Eq. 0- Thus multiplying Eq. jHoj i by B t = U _j Z7 // ; we find 

5ui m Bi [SXiUi - 5Rui - hi) 

= Bi (SRui+hi). 


( 60 ) 

□ 


B. Proof of Theorem [3] 

Proof. In order to prove Theorem^ we compute E[<5u ; (P) Su ^(P)] and E[5«; (P) 6uJ n (P)] based on the expression 
of ()U/ (P) which we found in Theorem [2] Then, results from |2j are used to simplify the computed expression. For 
convenience, we drop the dependency on P from Ui(P) and hfP) throughout the proof. 

Using the result from Theorem^ E[<5iq(P) 5u^ l {P)\ and E[<5{q(P) Su^ n (P)] are written as 

E[<MP) 6uZ(P)] « BfiKSRm + h t )(u^SR + h*)\B* 

= E {BiSRuiu^SRB^] + Bihih^B^ 


and 

E[5«z(P) Sul(P)] » BtEKSRm + hi){u T m 5R + h T m )\B T m 
= E [BtSRutulSRB*] + BAhlB T m . 


Using the following results from Theorem 9.2.2 in 


HE 


E [BtSRumZSRBZ] 


A; A i H 

i^l 


1 See also the proof of the Theorem 9.2.2 ^2 p. 454], 


(62) 
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and 


E [BiSRuiu^dRB^] 


in Eq. (f6T| and Eq. ([62)> proves the theorem. 


A; Am UjU ^ 

N (Xl-Xmf 
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